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Abstract 

In this article we present an artifact due to the conventional use of periodic boundary 
condition (PBC) which produces misleading results on charge transport behavior 
at low field regime. When simulated using PBC, artifact produces a negative field 
dependence of mobility (NFDM) at low field regime even when the positional dis- 
order is neglected. The artifact is due to the extra bias the carrier gains due to 
the neglected hops at boundaries which it encounters perpendicular to the field di- 
rection upon implementing the PBC. An alternative approach for simulating the 
charge transport using PBC without the observed artifact is also proposed. 

Key words: Molecularly doped polymers, Disordered system, Electronic transport, 
Negative field dependence of mobility, Monte Carlo simulation, Periodic boundary 
condition 

PACS: 72.80.Le, 72.80.Lg, 71.55.Jv, 73.61.Ph, 68.55.-a 



* Corresponding author. Tel.: +91 7312488361; Fax: -^91 7312488300 
Email address: rajmho@gmail.com (S. Raj Mohan). 



8 September 2009 



1 Introduction 



Gaussian disorder model (GDM) [1-4] is a widely accepted model that could 
explain most of the charge transport behavior observed in disordered organic 
systems. According to GDM the charge transport in disordered organic materi- 
als occur by hopping among transport sites that are subjected to energetic and 
positional disorder [1-4] . Since the model assumes Gaussian density of states 
[1-4] a complete analytical solution of the hopping transport is therefore dif- 
ficult, especially in 3-D. Hence the predictions of the GDM are made on the 
basis of Monte Carlo simulation of hopping charge transport [1-4]. Monte Carlo 
simulation is considered as an idealized experiment with which one can study 
the charge transport in disorder system as function of several parameters. Gen- 
erally charge transport is simulated for a sample length of few microns. This 
is to make sure that the carrier has attained a dynamic equilibrium during its 
transit and also to have a better comparison with the experiment, like Time of 
flight (TOF), which is generally performed on micron size thick samples [1-4]. 
In order to simulate the charge transport for a sample length of several microns 
periodic boundary condition (PBC) is frequently employed [la, 5-10]. The ad- 
vantage of using the PBC is that the simulation can be performed on a sample 
length of several microns using an array of smaller size. If PBC is not employed 
then an array of bigger sizes, in the all the three Cartesian directions, is re- 
quired which demands large computational resources. Even though the PBC 
have been used frequently in Monte Carlo simulation of charge transport in 
disordered organic systems a clear explanation of the method of implementing 
PBC in these simulation is scarce. Earlier literature [la, relevant references 
therein] reports that PBC have been implemented in the preliminary form 
as explained below. Implementing PBC in the preliminary form may lead to 
serious artifact which is the subject of this paper. 

The preliminary form of PBC is a set of boundary conditions that are used 
to simulate the properties of bulk system by simulating a part of it [11]. In 
principle PBC generates an inflnitely large system with help of a smaller array, 
that represent only a part of the bulk system, with the assumption that the 
small array will replicate periodically in all the three directions to form the 
bulk system. When we implement the PBC along one Cartesian direction, 
we take the carrier to the first plane of the array (lattice) in that direction 
when carrier reaches the final plane of the array in same direction. In this 
process the carrier's energy and other Cartesian coordinates remain same as 
at the boundary. When simulating the charge transport for several microns 
the carrier will encounter several such boundaries. Since the carrier is taken to 
the first plane when it reaches the final plane, in this process some of the hops 
that carrier may make in the absence of such a boundary get neglected, see 
Fig. 1(a). The neglected hops help the carrier to cover the required thickness 
in short time, which effectively provides a bias in addition to the applied field. 
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Number of such neglected hops at boundary can be significant particularly at 
smaller field strengths, higher temperature and higher energetic disorder. It is 
possible that these neglected hops can lead to an artifact which can seriously 
affect the charge transport simulation studies using PBC. This artifact, due 
to the neglected hops, can be significant especially at low electric field regime 
and is very critical as far as the investigations related to the mechanism of 
charge transport and the operation of devices like organic photovoltaic devices 
is concerned. 

In this article, we present our Monte Carlo simulation studies of charge trans- 
port on the basis of GDM [1-4] which unambiguously proves the artifact due 
to conventional use of PBC. At low field regime, field dependence of mobility 
simulated using PBC with zero positional disorder shows negative field de- 
pendence of mobility (NFDM) for all values of energetic disorder while a clear 
saturation of mobility with field was observed when simulated for the same 
inputs without using PBC. The observed NFDM in the absence of positional 
disorder contradicts with GDM [1-4,12] which asserts that the NFDM at low 
field regime can occur only in the presence of high positional disorder. Thus, 
the observed NFDM in the absence of any positional disorder is an artifact 
due to the conventional method of implementing the PBC. The origin of the 
artifact is attributed to the neglect of hops that the carrier may make in the 
absence of boundary which is created upon implementing the PBC. These 
neglected hops gives an extra bias to the carrier and there by enhancing the 
mobility. In concise, this study not only highlights an artifact, which can mis- 
lead interpretation and modeling of charge transport, but also cautions on the 
use of Monte Carlo simulation along with PBC for investigating the origin of 
NFDM / charge transport at low field regime. On the basis of the origin of 
artifact an alternative simulation approach for implementing the PBC which 
is free from the observed artifact is also proposed. 



2 Details of Monte Carlo simulation 

The Monte Carlo simulation is based on the commonly used algorithm re- 
ported by Schonherr et al [13]. A 3D array was considered as the lattice with 
size 70x70 along x and y direction. Along z direction, the direction of the 
applied field, various sizes were used to implement different simulation ap- 
proaches adopted for covering the required sample thickness along this direc- 
tion. Three different approaches are: 

Case 1 : In this lattice of size 70x70x70 along x, y and z direction 

was used for simulation (Fig. 1(b)). PBC was implemented along all the three 
directions. PBC along the applied field direction was implemented by taking 
the carrier to the plane z—1 when the carrier reaches 2;=70*'* plane, keeping 
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the other coordinates and the carrier energy same as at the boundary. PBC 
was also implemented along x and y direction in a similar manner. The dis- 
advantage of implementing the PBC in this manner is that some of the hops 
that the carrier may make in the absence of boundaries, which it encounters 
in the process of transit, along x, y and z directions are neglected. Simulations 
were also performed in a similar way with bigger lattice size along z direction 
(70x70x300, 70x70x1600, etc.). Size of the lattices was always chosen after 
considering the available computational resources. 

Case 2: In this case the simulation was performed without using PBC along 
z direction (Fig. 1(a)). Here the size of lattice along z direction was taken to 
be the sample length, which requires array of bigger size. However, PBC along 
X and y direction was implemented as explained above. In this case, carrier 
does not see any boundary along z direction and hence neglected hops at such 
boundaries in case 1 are taken into consideration. This simulation approach 
requires large computational resources. 

Case 3: In this case we show an effective way to implement the PBC for 
simulating the charge transport (Fig. 1(c)). In order to show an effective way 
to use PBC, a array of size 70x70x150 along x, y and z direction was used. 
Justification for the size of array is given later. Carrier is first injected into 
z=l plane and is allowed to move in the direction of applied electric field. The 
carrier was taken into z=70*^ plane when it reaches z=lAO^^ plane, keeping 
the x and y coordinates same. The energy of the carrier was kept same as at 
2=140*'* plane. PBC along x and y direction was implemented as mentioned 
before. Compared to case 1, this approach considers all those hops that are 
neglected at the boundaries perpendicular to z direction. Moreover, it requires 
less computational resources compared to case 2. In this case, when the PBC 
is implemented the carrier can perform all the hops around ^=70*^^ plane and 
hence the artifact due to neglected hops can be removed. 

Due to the PBC along x and y directions, in all the three cases, some hops may 
be still neglected when the carrier encounters a boundary perpendicular to x 
I y directions. In order to understand the influence of those neglected hops on 
charge transport simulation was performed, for all the three cases, by taking 
the carrier to the middle of the same plane (same z coordinate) when the 
carrier encounters a boundary either perpendicular to x / y direction. Carrier 
energy was kept same as that at the boundary. Here the carrier is allowed to 
make all the neglected hops around the middle of lattice. We found that the 
hops that may be neglected at the boundaries perpendicular to x / y direction 
upon implementing the PBC have negligible effect on charge transport. Data 
presented in this article, for all the three cases, was simulated using PBC along 
X / y direction. 

In this study simulation was always performed for a sample length of 4//m 
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along field direction. The lattice constant a — GA was used for the whole 
set of simulation [1] . The site energies of lattice were taken randomly from a 
Gaussian distribution with a known standard deviation (cr). Through out the 
simulation the positional disorder was neglected (E=0). This is to avoid the 
huge computer time required for simulating the charge transport with non-zero 
positional disorder. More over, wc believe that the outcome of this study can 
be unambiguously conveyed and justified even with present simulation with 
zero positional disorder. Simulation was always performed on this energetically 
disordered lattice with the assumption that the carrier hops among the lattice 
sites following Miller- Abrahams equation [14]. Throughout the simulation we 
took 270 = 10 [1,13]. Transit time of a carrier was calculated by adding 
all the hopping times and averaging over few hundred carriers. The mobility 
was calculated using drift mobility equation. The electric field range (~> 10^ 
V/cm) over which the simulations were performed is well above the field range 
(~10^-10^V/cm) over which the pure diffusion dominates the charge transport 
[12,15]. Hence the use of drift mobility equation over the field range of this 
study is appropriate. Simulation was preformed for various values of energetic 
disorder {a) and electric field strengths. All the data presented in this article 
was simulated with uncorrelated site energies. 



3 Results and Discussions 



Fig. 2 compares the simulated field dependence of mobility, at T=300K, for 
various values of energetic disorder, for easel and case 2. For both easel 
and case 2. field dependence of mobility for all the values of energetic disorder 
under study is similar except at low electric field strengths (~<3.6xl0^ V/cm). 
In easel, the field dependence of mobility at lower electric field strengths 
first decreases with increase of electric field and reaches a minimum value of 
mobility (represented by arrows in Fig. 2) before it shows positive dependence 
in a XogiiVs-E^^"^ fashion as predicted by GDM [1-4]. In case 2, for all the 
values of energetic disorder under study, the field dependence of mobility at 
lower electric field strengths show a clear saturation of mobility before it show 
positive dependence in a log fiVs.E^^'^ fashion as predicted by GDM [1-4]. Here 
the important point to highlight is that in easel we have observed the NFDM, 
at lower electric field strengths, even in the absence of positional disorder. 
This is against the GDM which strongly asserts that NFDM can occur only in 
the presence of high positional disorder [1-4,12]. When the energetic disorder 
decreases the NFDM at lower electric field strengths become remarkable. The 
strength of NFDM is assigned as the difference between the mobility value 
for the lowest electric field strength under study and the observed minima of 
the mobility at low field regime (shown by arrows in Fig. 2). The strength 
of NFDM increases with decrease of disorder as shown in Fig. 3. Here we 
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concentrate only on the origin of the difference in field dependence of mobility 
at lower electric field strengths between easel and case 2. Simulation similar 
to case 1 reported by H. Basslcr [la] shows a saturation of mobility at low 
field regime. Unfortunately, in the reported data presents very few number 
of data points at low field regime which is not all sufficient to draw a clear 
conclusion about the nature of field dependence of mobility. Thus the authors 
might have missed the very small NFDM present at high value of energetic 
disorder (~0.1eV) which is generally used for simulation. 

In both case 1 and case 2 the simulation was performed for same input pa- 
rameters. The only different is that in case 1 the required sample length (4/im) 
is covered with use of PBC and in case 2 without using PBC. So the observed 
difference in field dependence of mobility at low field regime between case 1 
and case 2 can be an artifact of PBC. When PBC is implemented as in case 
1, the carrier is taken to the first plane {z—1 or x—1 or y—1, depending on the 
direction) when it reaches the final plane {z=70 or x=70 or y=70). In that 
process some hops that the carrier may make in the absence of a boundary 
can be neglected. In the absence of a boundary the carrier may wander more 
before proceeding further in the field direction (Fig. la). So in case 1, some 
hops are neglected while implementing the PBC. Therefore, the number of 
hops made by the carrier in case 2 must be higher than in case 1. Compared 
to high electric field strength, carrier wanders more in the lattice at low elec- 
tric field strength during its transit along the field direction. Thus the number 
of hops neglected at the boundaries upon implementing the PBC is expected 
to be high at low field and vice versa. This predicts a remarkable difference 
between the total number of hops made by the carrier in case 1 and 2 and 
this difference is expected to be high at lower electric field strength. This 
difference should gradually decreases with increase of electric field strength 
and become negligibly small at higher electric field strengths. Compared to 
low value of energetic disorder the carrier wanders more in the lattice at high 
value of energetic disorder during its transit. At low value of energetic disorder 
carrier is forced to move mostly in the direction of the applied field. Similar 
to the above explanation, for a constant field strength (for example the lowest 
field strength used in the study, 4x10"^ V/cm) the difference between the total 
number of hops, made by the carrier in case 1 and 2, should increase with 
the increase of energetic disorder. We have confirmed these conjectures in our 
simulation. 

Fig. 4 shows the variation of the total number of hops made by the carrier in 
case 1 and 2, at T=300K, with electric field for various values of energetic dis- 
order. Fig. 4 clearly shows that at low electric field strengths, for all the values 
of energetic disorder under study, the total number of hops made by the carrier 
in case 2 is higher than in case 1. This difference between the total number 
of hops for the two cases {case 1 and case 2) is higher for the lowest electric 
field strength used in the study (4x10^ V/cm). This difference decreases with 



6 



increase of electric field strength and become negligible at very high electric 
field strengths. It also shows that the difference in the total numbers of hops 
(hops), made by the carrier in two cases {case 1 and is higher for higher 
value of energetic disorder (inset of Fig. 4). Thus, Fig. 4 unambiguously estab- 
lishes the fact that upon implementing the PBC substantial number of hops 
that the carrier may make in the absence of boundary is neglected especially 
at low electric field strengths as well as at high value of energetic disorder. 

The neglect of hops, in case 1, upon implementing the PBC certainly influ- 
ences the charge transport especially the transit time. It is confirmed through 
simulation that in both cases the charge carriers have attained dynamic equi- 
librium [8,16] in a similar fashion while covering the required sample length. 
So the neglected hops do not affect the energy relaxation of carrier during the 
hopping in Gaussian DOS. Fig. 5 shows the field dependence of the difference 
between the carrier transit time (Ar) for case 2 and case 1 for various value 
of energetic disorder. For all values of energetic disorder under study, the max- 
imum difference in carrier transit time between case 2 and case 1 is observed 
for the lowest electric field strength. This difference gradually decreases with 
increase of electric field strength and become neghgibly small at higher elec- 
tric field strengths. In case 1, the neglect of hops effectively gives an extra 
bias for the carrier to move in the applied field direction. In principle, carrier 
covers the required sample length in less number of hops. Thus the transit 
time of the carrier in case 1 must be less than that of case 2. The maximum 
difference in transit time between case 1 and case 2 must occur for the case 
with maximum difference in the average number of hops. Hence the maximum 
difference in transit time must occur at low field strengths and this difference 
should gradually decreases with increase of electric field strength which ulti- 
mately become negligibly small at high electric field strengths (Fig. 5). For the 
same reason the difference in transit time between case 1 and case 2 must be 
higher for higher value of energetic disorder (Fig. 5) and this difference should 
decrease with decrease of energetic disorder. Having shown how the difference 
in transit time, for case 1 and case 2, varies as a function of electric field and 
energetic disorder, we would like to relate the difference in transit time to the 
difference in mobility, 

L 

Hi-H2 = — 

where L is the thickness of sample, E is the applied electric field, ^ = At 
is the difference in transit time between case 1 and case 2, T2 is transit time 
in case 2 and Ti is the transit time in case 1. Equation 1 suggests that the 
difference in mobility, for a constant electric field and thickness, depends on 
the value of ^ rather than on Ar alone. In Fig. 6 electric field is limited to the 
range where the value of Ar is significant. For any value of energetic disorder, 
^ decreases with increase of electric field strengths (Fig. 6(a)). This supports 



i 
ri 



i 

To, 



L 



(1) 



7 



the observed NFDM at low field regime for all values of energetic disorder 
under study. It is clear from Fig. 6(b) that the value of ^ also decreases with 
increase of energetic disorder. This suggests a strong NFDM for low value of 
energetic disorder and vice versa as shown in Fig. 3. 

The above presented data and its explanation clearly showed that an artifact 
can occur at low electric field regime upon implementing the PBC in conven- 
tional way. The observed artifact is reduced when simulated in a lattice with 
large size along the field direction. When the lattice size along the field direc- 
tion is large, carrier encounters fewer boundaries while covering the required 
sample length. This reduces the number of neglected hops and hence the ar- 
tifact (data not shown). This observation stresses the influence of neglected 
hops and supports the explanation of observed artifact in case 1. This method 
is inconvenient because the artifact cannot be completely removed and big- 
ger arrays require large computations resources as well. More over, the size of 
the array required depends on chosen sample length and energetic disorder. 
In case 3 a different method of simulation is adopted with which the PBC 
can be implemented without the observed artifact. Hence the simulation was 
performed on a lattice of size 70x70x150 along x, y and z direction. The car- 
rier was injected into the first plane {z—1) and allowed to move till 2; =140*'* 
plane. Once it reaches the z=14Q^^ plane the carrier was taken to z—lQ plane, 
keeping the energy of carrier and x and y co-ordinates the same as at z=140*'* 
plane. In case 3, the carrier is allowed to make all the required number of hops 
around z=70 plane when the PBC is used. Thus the hops that are neglected 
at the boundaries in case 1 are properly accounted for, there by eliminating 
the observed artifact. Therefore, the field dependence of mobility and the to- 
tal number of hops made by the carrier in case 3 is expected to be similar 
as in case 2 for all the values of energetic disorder. This is clearly shown in 
from Fig. 7 where the data in both cases {case 2 and case 3) super impose 
each other. Fig. 7 asserts the fact that NFDM observed in case 1 is due to 
extra bias gained by the carrier due to the neglected hops. Care should be 
taken while choosing the dimension of the array because carrier in any case 
should not reach the z=l plane from the plane it has taken (2:=70 in our case) 
after the use of PBC. If carrier move back and touches z=l plane, then some 
hops may be neglected and method become less effective which in turn may 
result in NFDM. Hence a sufficient buffer lattice must be provided for the 
wandering of the carrier around the plane where it has taken after the use of 
PBC. In this study, the optimum buffer size (70x70x70) is obtained by com- 
paring the field dependence of mobility for different buffer sizes with the field 
dependence of mobility for similar inputs using case 2 (data not shown). The 
presence of small buffer between z=140*'* and z=150*'* plane makes sure that 
the probability of jumping is always calculated for the same number of sites. 
Our simulation studies also showed that PBC applied along x and y direction 
has negligible effect on the field dependence of mobility (data not shown). 
This suggests that the observed artifact is due to the substantial number of 
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hops that arc neglected at the boundaries perpendicular to the field direction 
{z direction) upon implementing the PBC. The reason may be that the num- 
ber of boundaries carrier encounters perpendicular to x and y direction is less 
compared to the number of boundaries the carrier encounters perpendicular 
to z direction, which is inevitable while covering the required sample length. 



4 Conclusion 

This study exposes an artifact in conventional method of implementing the 
PBC in Monte Carlo simulation of charge transport in disordered organic sys- 
tems. This manifests in NFDM at low field regime for all values of energetic 
disorder under study even in the absence of positional disorder. Upon imple- 
menting the PBC a substantial number of hops that the carrier may make in 
the absence of boundary, which it encounters perpendicular to the field direc- 
tion, are neglected. Effectively, the carrier covers the required sample length 
in less number of hops resulting in a shorter transit time. Thus the carrier 
gains an extra bias along the field direction. Hence the origin of artifact is 
rationalized on the basis of this extra bias gained by the carrier to move in 
the applied field direction. The simualtions confirmed that the boundary con- 
dition applied in directions other than applied field directions have negligible 
effect on charge transport. An alternative approach for simulation with which 
one can implement the PBC without the observed artifact is also proposed. In 
cocise this study cautions the researchers who adopt periodic boundary condi- 
tion in Monte Carlo simulation for studying the charge transport mechanism 
especially at low field regime. 
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Figure Captions 



Figure 1. (a) Schematic diagram showing a typical hopping motion of car- 
rier inside the lattice. Shaded region shows the hops that are neglected upon 
implementing the PBC. (b) Schematic diagram showing case 1 (c) schematic 
diagram showing case 3. Schematic diagram of case 2 is same as shown in (a). 

Figure 2. Comparison of field dependence of mobility, at T=300K, for vari- 
ous values of disorder for (o) case 1 and (•) ca.se 2. Arrow shows the minima 
of mobility occurred at low field regime in case 1. Inset shows the magnified 
view of low field regime for the respective cases. 

Figure 3. Variation of the strength of negative field dependence of mobility, 
observed in case 1, with energetic disorder, at T=300K. Solid fine shows a 
guide to eye. 

Figure 4. Field dependence of the total number of hops made by the carrier 
in two cases {case 1 (solid line) and case 2 (dashed line)), for various values 
of energetic disorder, at T=300K. Inset shows the dependence of difference in 
total number of hops (hops), made by the carrier in two cases {case 1 and 
case 2), on energetic disorder at 4xlO'^V/cm and T=300K. 

Figure 5. Field dependence of the difference between the carrier transit time 
(At), for case 1 and case 2, for various values of energetic disorder at T=300K. 

Figure 6. Variation of ^ (a) with electric field strength for the values ener- 
getic disorder where reasonably strong NFDM is observed (b) with energetic 
disorder for the various value of electric field strength at low field regime. Solid 
line in (b) is a guide to eye. Averaging is carried over 500 carriers at T=300K. 

Figure 7. Comparison of the field dependence of mobility, at T=300K, simu- 
lated for (•) case 2 and (o) case 3. Inset shows the comparison of total number 
of hops made by carrier in (a) case 2 and (o) case 3. 
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Figure 1 
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Figure 3 
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Figure 4 
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